function [x,Q,diffQ,idx,residuals,yfit,numObs] = NLSfilter(y,model,x0,lambda0,tolFun,t)
%This function computes the NLS filter where optimization is done 
%by the Levenberg Marquardt algorhtm. 
display    = 0;
if model.factorSignResOn == 1
    % Transform x0
    x = [log(x0(1:model.nm));x0(model.nm+1:end,1)];
else
    x = x0;
end
nx         = size(x0,1);
epsilon    = 1D-7;
yfit_p_eps = zeros(size(y,1),nx);
yfit_m_eps = zeros(size(y,1),nx);
diffQ      = 1e6;
idx        = 0;
lambda     = lambda0;
updateJ    = 1;
yfit       = modelFit(model,x,t);
select     = ~isnan(y);
numObs     = sum(select,1);
residuals  = y - yfit;
Q          = sum(residuals(select,1).^2,1);
while diffQ > tolFun && idx < 100
    if updateJ == 1
        J = getJacobianNLS(model,x,t);
        %for i=1:nx
        %    x_p_eps = x;
        %    x_p_eps(i,1) = x_p_eps(i,1) + epsilon;
        %    yfit_p_eps(:,i) = modelFit(model,x_p_eps,t);
        %    
        %    x_m_eps = x;
        %    x_m_eps(i,1) = x_m_eps(i,1) - epsilon;
        %    yfit_m_eps(:,i) = modelFit(model,x_m_eps,t);
        %end
        %Jnum = (yfit_p_eps(select,:) - yfit_m_eps(select,:))/(2*epsilon);
        %max(abs(J(:)-Jnum(:)))
        %J = Jnum;
    end
    
    % Computing the new value of x
    %delta = (J'*J + lambda*diag(diag(J'*J)))\J'*(y - yfit);
    delta = (J'*J + lambda*eye(nx))\J'*(y(select,1) - yfit(select,1));
    
    % Evaluating fit at new x
    x_new         = x + delta;
    yfit_new      = modelFit(model,x_new,t);
    residuals_new = y - yfit_new;
    Q_new         = sum(residuals_new(select,1).^2,1);
    diffQ         = abs(Q_new - Q);
    if display == 1
        disp(['Q           = ', num2str(Q_new)])
        disp(['lambda      = ', num2str(lambda)])
        disp(['change in Q = ', num2str(Q_new - Q)]);
        disp(['x_new =', num2str(x_new')])
        disp(['delta =', num2str(delta')])
        
    end
    if Q_new < Q
        lambda = lambda/10;
        x = x_new;
        residuals = residuals_new;
        yfit = yfit_new;
        Q = Q_new;
        updateJ = 1;
    else
        updateJ = 0;
        lambda = lambda*10;
    end
        
    % Update index
    idx   = idx +1;
end

% Untransform x;
if model.factorSignResOn == 1
    x = [exp(x(1:model.nm,1));x(model.nm+1:end,1)];
end

    
end
